% return a quaternion from degree and rotate axis
function q = quaternion(etha, n)
  etha /= 2;
  q = [cosd(etha), n(1)*sind(etha), n(2)*sind(etha), n(3)*sind(etha)];
endfunction